EDA: Cactus Flowering Season Exploratory Data Analysis

Here, we will explore the data collected during the Spring 2019 flowering field season at the Granite Mountains desert research center.

Some key datasets: * Bird visitation: the number of bird visitations to a focal cactus individual, with manipulated showiness of flowers (both real and fake) on a cactus at zero, fifteen, and thirty flowers per cactus. First dataset is only focal 1 hour observations. * Cactus allocation: a measurement of reproductive and non reproductive structures of 100 cactus individuals. This dataset will be paired with data on fruit production of the same individuals. * Side Wide: entire-site measurements of bird diversity, mesohabitats, and behaviors collected by line transect walks * Cactus architecture temperature: records of micro-habitat temperature related to cacti by measuring the top of the cactus’s canopy, within the cactus’s canopy, and the cactus’s understory.

#libraries to include
library(ggplot2)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(kableExtra)
library(ggpubr)

## Loading required package: magrittr

library(tidyverse)

## -- Attaching packages ----------------------------------------------------------------- tidyverse 1.2.1 --

## v tibble  1.4.2     v purrr   0.2.5
## v tidyr   0.8.1     v stringr 1.3.1
## v readr   1.1.1     v forcats 0.3.0

## -- Conflicts -------------------------------------------------------------------- tidyverse_conflicts() --
## x tidyr::extract()   masks magrittr::extract()
## x dplyr::filter()    masks stats::filter()
## x dplyr::lag()       masks stats::lag()
## x purrr::set_names() masks magrittr::set_names()

library(reshape)

## 
## Attaching package: 'reshape'

## The following objects are masked from 'package:tidyr':
## 
##     expand, smiths

## The following object is masked from 'package:dplyr':
## 
##     rename

library(reshape2)

## 
## Attaching package: 'reshape2'

## The following objects are masked from 'package:reshape':
## 
##     colsplit, melt, recast

## The following object is masked from 'package:tidyr':
## 
##     smiths

library(ggmap)

## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.

## Please cite ggmap if you use it! See citation("ggmap") for details.

## 
## Attaching package: 'ggmap'

## The following object is masked from 'package:magrittr':
## 
##     inset

library(DT)

Site Wide Exploration

Here, we will primarily examine the behavior, the biodiversity, and mesohabitat use of/by all birds present in the entire site (that is, excluded birds that were flying over, like Turkey Vultures soaring over). This dataset functions as a large control for our focal-cactus bird visitation observations.

#bring in dataset of interest
site <- read.csv("~/Masters/Flower_CactusBirdInteractionR/data/line_transects/line_transects.csv")
head(site)

#summary stats
unique(site$species)

##  [1] Mourning Dove                         
##  [2] Canyon Wren                           
##  [3] Phainopepla                           
##  [4] House Finch                           
##  [5] Unknown                               
##  [6] Violet.green Swallow                  
##  [7] Hummingbird                           
##  [8] Northern Mockingbird                  
##  [9] Black.throated Sparrow                
## [10] Blue.gray Gnatcatcher                 
## [11] Black.headed Grosbeak                 
## [12] Ash.throated Flycatcher               
## [13] Wilson's Warbler                      
## [14] Common Raven                          
## [15] Black.tailed Gnatcatcher              
## [16] Rock Wren                             
## [17] Costa's Hummingbird                   
## [18] Say's Phoebe                          
## [19] House finch                           
## [20] Green.tailed Towhee                   
## [21] Crissal's Thrasher                    
## [22] Hammond's Flycatcher                  
## [23] Black.chinned hummingbird             
## [24] Anna's Hummingbird                    
## [25] White.throated Swift                  
## [26] MacGillivray's Warbler                
## [27] Gambel's Quail                        
## [28] Verdin                                
## [29] Cactus Wren                           
## [30] Western Wood.Pewee                    
## [31] Great Horned Owl                      
## [32] Unknown Passerine                     
## [33] Gray Flycatcher                       
## [34] Pacific.slope Flycatcher              
## [35] Townsend's Warbler                    
## [36] Warbling Vireo                        
## [37] Hooded Oriole                         
## [38] Nuttall's/Ladderback Woodpecker Hybrid
## [39] Red.tailed Hawk                       
## [40] Western Tanager                       
## [41] Nashville Warbler                     
## 41 Levels: Anna's Hummingbird ... Wilson's Warbler

#41 unique species
unique(site$behavior)

##  [1] Singing                 Flying                 
##  [3] Standing on Ground      Flying Between Plants  
##  [5] Sitting on Plant        Flying on Plant        
##  [7] Foraging                Calling                
##  [9] Sparring                Carrying Nest Materials
## [11] Copulating              Walking                
## [13] Distraction Display     Preening               
## 14 Levels: Calling Carrying Nest Materials ... Walking

#14 unique behaviors
unique(site$mesohabitat)

##  [1] no visual         air               rock             
##  [4] cactus canopy     shrub             cactus understory
##  [7] gravel            grass             annuals          
## [10] wash             
## 10 Levels: air annuals cactus canopy cactus understory grass ... wash

#10 unique mesohabitats


#Alternate dataset that only include bird visuals
site_visuals <- na.omit(site)
head(site_visuals)

#summary stats
unique(site_visuals$species)

##  [1] Phainopepla                           
##  [2] Mourning Dove                         
##  [3] House Finch                           
##  [4] Unknown                               
##  [5] Violet.green Swallow                  
##  [6] Hummingbird                           
##  [7] Northern Mockingbird                  
##  [8] Black.throated Sparrow                
##  [9] Blue.gray Gnatcatcher                 
## [10] Black.headed Grosbeak                 
## [11] Ash.throated Flycatcher               
## [12] Wilson's Warbler                      
## [13] Black.tailed Gnatcatcher              
## [14] Rock Wren                             
## [15] Costa's Hummingbird                   
## [16] House finch                           
## [17] Green.tailed Towhee                   
## [18] Crissal's Thrasher                    
## [19] Hammond's Flycatcher                  
## [20] Black.chinned hummingbird             
## [21] Say's Phoebe                          
## [22] Anna's Hummingbird                    
## [23] White.throated Swift                  
## [24] MacGillivray's Warbler                
## [25] Gambel's Quail                        
## [26] Verdin                                
## [27] Unknown Passerine                     
## [28] Gray Flycatcher                       
## [29] Pacific.slope Flycatcher              
## [30] Warbling Vireo                        
## [31] Townsend's Warbler                    
## [32] Cactus Wren                           
## [33] Hooded Oriole                         
## [34] Nuttall's/Ladderback Woodpecker Hybrid
## [35] Western Wood.Pewee                    
## [36] Nashville Warbler                     
## 41 Levels: Anna's Hummingbird ... Wilson's Warbler

#41 unique species
unique(site_visuals$behavior)

##  [1] Flying                  Standing on Ground     
##  [3] Flying Between Plants   Singing                
##  [5] Sitting on Plant        Flying on Plant        
##  [7] Foraging                Calling                
##  [9] Sparring                Carrying Nest Materials
## [11] Copulating              Walking                
## [13] Distraction Display     Preening               
## 14 Levels: Calling Carrying Nest Materials ... Walking

#14 unique behaviors
unique(site_visuals$mesohabitat)

##  [1] air               rock              cactus canopy    
##  [4] shrub             cactus understory gravel           
##  [7] grass             no visual         annuals          
## [10] wash             
## 10 Levels: air annuals cactus canopy cactus understory grass ... wash

#10 unique mesohabitats
#counts are the same for visuals versus no visuals

#map of all bird sitings
deep.sunset <- get_map(location = c(lon = -115.663, lat = 34.7805), zoom = 17, color="bw")

## Source : https://maps.googleapis.com/maps/api/staticmap?center=34.7805,-115.663&zoom=17&size=640x640&scale=2&maptype=terrain&language=en-EN&key=xxx-6x1tmeuVFVYhc8MMkmuUuKOALLNc

site.wide.map <- ggmap(deep.sunset)
site.wide.map <- site.wide.map +
  geom_point(data=site, aes(x=lon, y=lat, colour = family, group = family), alpha = 3/10, size =4) +
  labs(title = "Site Wide Species Observations", x = "longitude", y = "latitude", color = "Family")
site.wide.map

#histogram of bird behaviors
behav <- ggplot(site, aes(x=behavior)) + geom_bar() + theme(axis.text.x = element_text(angle = 60, hjust = 1))
#lump into smaller groups 
bbehav <- ggplot(site, aes(x=broad_behavior)) + geom_bar() + theme(axis.text.x = element_text(angle = 60, hjust =1))

#histogram of bird mesohabitats
meso <- ggplot(site, aes(x=mesohabitat)) + geom_bar() + theme(axis.text.x = element_text(angle = 60, hjust = 1))

#histogram of bird species
species <- ggplot(site, aes(x=species)) + geom_bar() + theme(axis.text.x = element_text(angle = 60, hjust = 1))
#lump into broader taxa and functional groups
order <- ggplot(site, aes(x=order)) + geom_bar() + theme(axis.text.x = element_text(angle = 60, hjust = 1))
family <- ggplot(site, aes(x=family)) + geom_bar() + theme(axis.text.x = element_text(angle = 60, hjust = 1))

#histogram of bird distance from transect
dist <- ggplot(site, aes(x=distance)) + geom_bar() 

#histogram of bird behaviors
behavv <- ggplot(site_visuals, aes(x=behavior)) + geom_bar() + theme(axis.text.x = element_text(angle = 60, hjust = 1))
#lump into smaller groups 
bbehavv <- ggplot(site_visuals, aes(x=broad_behavior)) + geom_bar() + theme(axis.text.x = element_text(angle = 60, hjust =1))

#histogram of bird mesohabitats
mesov <- ggplot(site_visuals, aes(x=mesohabitat)) + geom_bar() + theme(axis.text.x = element_text(angle = 60, hjust = 1))

#histogram of bird species
speciesv <- ggplot(site_visuals, aes(x=species)) + geom_bar() + theme(axis.text.x = element_text(angle = 60, hjust = 1))
#lump into broader taxa and functional groups
orderv <- ggplot(site_visuals, aes(x=order)) + geom_bar() + theme(axis.text.x = element_text(angle = 60, hjust = 1))
familyv <- ggplot(site_visuals, aes(x=family)) + geom_bar() + theme(axis.text.x = element_text(angle = 60, hjust = 1))

#histogram of bird distance from transect
distv <- ggplot(site_visuals, aes(x=distance)) + geom_bar() 

#put together on one plot for side-by-side comparison

#taxanomic grid
tax <- ggarrange(species, speciesv, family, familyv, order, orderv,
                    labels = c("All", "Visuals", "All", "Visuals", "All", "Visuals"),
                    ncol = 2, nrow = 3)
tax

#difficult to view in grid arrangement
species

speciesv

family

familyv

order

orderv

#behavior grid
bhv <- ggarrange(behav, behavv, bbehav, bbehavv,
                   labels = c("All","Visuals", "All", "Visuals"),
                   ncol = 2, nrow = 2)
bhv

#mesohabitat grid
mesohabitat <- ggarrange(meso, mesov,
                         labels = c("All", "Visuals"),
                         ncol = 1, nrow =2)
mesohabitat

#distance grid
distance <- ggarrange(dist, distv,
                      labels = c("All", "Visuals"),
                      ncol = 1, nrow =2)
distance

#Behavior stats: which behaviors are most common?
#get counts for each behavior daily
behav <- (data.frame(table(site$behavior, site$date, site$walk)))

shapiro.test(behav$Freq)

## 
##  Shapiro-Wilk normality test
## 
## data:  behav$Freq
## W = 0.14823, p-value < 2.2e-16

#it's normally distributed

#run ANOVA 
bhv <- aov(Freq ~ Var1, behav)
summary(bhv)

##               Df Sum Sq Mean Sq F value Pr(>F)    
## Var1          13    453   34.82   16.08 <2e-16 ***
## Residuals   3626   7852    2.17                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#p-value < 2e-16, very significant (as expected)

#Now see what stood out and in what direction
#Tukey HSD
tuk.behav <- TukeyHSD(bhv, "Var1", ordered = TRUE)

#Make viewable
tuktab.behav <- data.frame(tuk.behav$Var1)
tuktab.behav <- tuktab.behav %>% rownames_to_column("behavior") 
tuktab.behav <- separate(tuktab.behav, behavior, into = c("behav1", "behav2"), sep = "-")
tuktab.behav$diff <- NULL
tuktab.behav$lwr <- NULL 
tuktab.behav$upr<- NULL 
# Make into "wide" format using Var2
tuktab.behav <- spread(data = tuktab.behav, key = 'behav2', value = 'p.adj', fill = '')
tuktab.behav

#Would love to figure out how to color cells conditionally....

#Do it for visuals only
behav.v <- (data.frame(table(site_visuals$behavior, site_visuals$date, site_visuals$walk)))

shapiro.test(behav.v$Freq)

## 
##  Shapiro-Wilk normality test
## 
## data:  behav.v$Freq
## W = 0.16204, p-value < 2.2e-16

#it's normal

bhvv <- aov(Freq ~ Var1, behav.v)
summary(bhvv)

##               Df Sum Sq Mean Sq F value   Pr(>F)    
## Var1          13   63.9   4.918   8.079 3.18e-16 ***
## Residuals   3626 2207.0   0.609                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#significant model

#Now see what stood out and in what direction
#Tukey HSD
tukbehav.v <- TukeyHSD(bhvv, "Var1", ordered = TRUE)

#Make viewable
tuktab.behav.v <- data.frame(tukbehav.v$Var1)
tuktab.behav.v <- tuktab.behav.v %>% rownames_to_column("behavior") 
tuktab.behav.v <- separate(tuktab.behav.v, behavior, into = c("behav1", "behav2"), sep = "-")
tuktab.behav.v$diff <- NULL
tuktab.behav.v$lwr <- NULL 
tuktab.behav.v$upr<- NULL 
# Make into "wide" format using Var2
tuktab.behav.v <- spread(data = tuktab.behav.v, key = 'behav2', value = 'p.adj', fill = '')
tuktab.behav.v

#Do we still see significant differences when we lump similar behaviors together?
#Do it for visuals only
broad.behav <- (data.frame(table(site$broad_behavior, site$date, site$walk)))

shapiro.test(broad.behav$Freq)

## 
##  Shapiro-Wilk normality test
## 
## data:  broad.behav$Freq
## W = 0.17813, p-value < 2.2e-16

#it's normal

bb <- aov(Freq ~ Var1, broad.behav)
summary(bb)

##               Df Sum Sq Mean Sq F value Pr(>F)    
## Var1           5    793  158.61   17.58 <2e-16 ***
## Residuals   1554  14018    9.02                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#significant model

#Now see what stood out and in what direction
#Tukey HSD
tuk1.bb <- TukeyHSD(bb, "Var1", ordered = TRUE)

#Make viewable
tuktab.behav.bb <- data.frame(tuk1.bb$Var1)
tuktab.behav.bb <- tuktab.behav.bb %>% rownames_to_column("behavior") 
tuktab.behav.bb <- separate(tuktab.behav.bb, behavior, into = c("behav1", "behav2"), sep = "-")
tuktab.behav.bb$diff <- NULL
tuktab.behav.bb$lwr <- NULL 
tuktab.behav.bb$upr<- NULL 
# Make into "wide" format using Var2
tuktab.behav.bb <- spread(data = tuktab.behav.bb, key = 'behav2', value = 'p.adj', fill = '')
tuktab.behav.bb

#How about broad behaviors but only visuals?
broad.behav.v <- (data.frame(table(site_visuals$broad_behavior, site_visuals$date, site_visuals$walk)))
shapiro.test(broad.behav.v$Freq)

## 
##  Shapiro-Wilk normality test
## 
## data:  broad.behav.v$Freq
## W = 0.20657, p-value < 2.2e-16

#it's normal

bbv <- aov(Freq ~ Var1, broad.behav.v)
summary(bbv)

##               Df Sum Sq Mean Sq F value   Pr(>F)    
## Var1           5    107  21.423   9.226 1.15e-08 ***
## Residuals   1554   3608   2.322                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#significant model

#Now see what stood out and in what direction
#Tukey HSD
tuk1.bb.v <- TukeyHSD(bbv, "Var1", ordered = TRUE)
#Make viewable
tuktab.bb.v <- data.frame(tuk1.bb.v$Var1)
tuktab.bb.v <- tuktab.bb.v %>% rownames_to_column("behavior") 
tuktab.bb.v <- separate(tuktab.bb.v, behavior, into = c("behav1", "behav2"), sep = "-")
tuktab.bb.v$diff <- NULL
tuktab.bb.v$lwr <- NULL 
tuktab.bb.v$upr<- NULL 
# Make into "wide" format using Var2
tuktab.bb.v <- spread(data = tuktab.bb.v, key = 'behav2', value = 'p.adj', fill = '')
tuktab.bb.v

#Difference in Species?
taxa.ssp <- (data.frame(table(site$species, site$date, site$walk)))
qqnorm(taxa.ssp$Freq)

#it's normal? Maybe not... cannot run shapiro due to too large sample size

sspa <- aov(Freq ~ Var1, taxa.ssp)
summary(sspa)

##                Df Sum Sq Mean Sq F value Pr(>F)    
## Var1           40    309   7.735   15.84 <2e-16 ***
## Residuals   10619   5187   0.488                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#significant model

#Now see what stood out and in what direction
#Tukey HSD
tuk2.ssp <- TukeyHSD(sspa, "Var1", ordered = TRUE)
#Make viewable
tuktab.ssp <- data.frame(tuk2.ssp$Var1)
tuktab.ssp <- tuktab.ssp %>% rownames_to_column("species") 
tuktab.ssp <- separate(tuktab.ssp, species, into = c("ssp1", "ssp2"), sep = "-")
tuktab.ssp$diff <- NULL
tuktab.ssp$lwr <- NULL 
tuktab.ssp$upr<- NULL 
# Make into "wide" format using Var2
tuktab.ssp <- spread(data = tuktab.ssp, key = 'ssp2', value = 'p.adj', fill = '')
#Kind of impossible to read, even in simplified format
tuktab.ssp

#Difference in Families?
taxa.fam <- (data.frame(table(site$family, site$date, site$walk)))
qqnorm(taxa.fam$Freq)

#it's normal? Maybe not... cannot run shapiro due to too large sample size

fama <- aov(Freq ~ Var1, taxa.fam)
summary(fama)

##               Df Sum Sq Mean Sq F value Pr(>F)    
## Var1          21    308  14.650   14.48 <2e-16 ***
## Residuals   5698   5764   1.012                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#significant model

#Now see what stood out and in what direction
#Tukey HSD
tuk2.fam <- TukeyHSD(fama, "Var1", ordered = TRUE)
#Make viewable
tuktab.fam <- data.frame(tuk2.fam$Var1)
tuktab.fam <- tuktab.fam %>% rownames_to_column("taxa") 
tuktab.fam <- separate(tuktab.fam, taxa, into = c("taxa1", "taxa2"), sep = "-")
tuktab.fam$diff <- NULL
tuktab.fam$lwr <- NULL 
tuktab.fam$upr<- NULL 
# Make into "wide" format using Var2
tuktab.fam <- spread(data = tuktab.fam, key = 'taxa2', value = 'p.adj', fill = '')
tuktab.fam

#Difference in Orders?
taxa.ord <- (data.frame(table(site$order, site$date, site$walk)))
shapiro.test(taxa.ord$Freq)

## 
##  Shapiro-Wilk normality test
## 
## data:  taxa.ord$Freq
## W = 0.10687, p-value < 2.2e-16

qqnorm(taxa.ord$Freq)

#it's normal

orda <- aov(Freq ~ Var1, taxa.ord)
summary(orda)

##               Df Sum Sq Mean Sq F value Pr(>F)    
## Var1           7   1694  241.99   23.98 <2e-16 ***
## Residuals   2072  20914   10.09                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#significant model

#Now see what stood out and in what direction
#Tukey HSD
tuk2.ord <- TukeyHSD(orda, "Var1", ordered = TRUE)
#Make viewable
tuktab.ord <- data.frame(tuk2.ord$Var1)
tuktab.ord <- tuktab.ord %>% rownames_to_column("order") 
tuktab.ord <- separate(tuktab.ord, order, into = c("order1", "order2"), sep = "-")
tuktab.ord$diff <- NULL
tuktab.ord$lwr <- NULL 
tuktab.ord$upr<- NULL 
# Make into "wide" format using Var2
tuktab.ord <- spread(data = tuktab.ord, key = 'order2', value = 'p.adj', fill = '')
#More passerines than anything else, for sure
tuktab.ord

#Let's try the visuals-only versions
#Species level
taxa.ssp.v <- (data.frame(table(site_visuals$species, site_visuals$date, site_visuals$walk)))
qqnorm(taxa.ssp.v$Freq)

#it's normal? Maybe not... cannot run shapiro due to too large sample size

sspav <- aov(Freq ~ Var1, taxa.ssp.v)
summary(sspav)

##                Df Sum Sq Mean Sq F value Pr(>F)    
## Var1           40   84.2  2.1049   10.94 <2e-16 ***
## Residuals   10619 2042.9  0.1924                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#significant model

#Now see what stood out and in what direction
#Tukey HSD
tuk2.ssp.v <- TukeyHSD(sspav, "Var1", ordered = TRUE)
#Make viewable
tuktab.ssp.v <- data.frame(tuk2.ssp.v$Var1)
tuktab.ssp.v <- tuktab.ssp.v %>% rownames_to_column("species") 
tuktab.ssp.v <- separate(tuktab.ssp.v, species, into = c("ssp1", "ssp2"), sep = "-")
tuktab.ssp.v$diff <- NULL
tuktab.ssp.v$lwr <- NULL 
tuktab.ssp.v$upr<- NULL 
# Make into "wide" format using Var2
tuktab.ssp.v <- spread(data = tuktab.ssp.v, key = 'ssp2', value = 'p.adj', fill = '')
tuktab.ssp.v

#Family visuals only
taxa.fam.v <- (data.frame(table(site_visuals$family, site_visuals$date, site_visuals$walk)))
qqnorm(taxa.fam.v$Freq)

#it's normal? Maybe not... cannot run shapiro due to too large sample size

famav <- aov(Freq ~ Var1, taxa.fam.v)
summary(famav)

##               Df Sum Sq Mean Sq F value Pr(>F)    
## Var1          21   85.5   4.071   10.27 <2e-16 ***
## Residuals   5698 2258.7   0.396                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#significant model

#Now see what stood out and in what direction
#Tukey HSD
tuk2.fam.v <- TukeyHSD(famav, "Var1", ordered = TRUE)
#Make viewable
tuktab.fam.v <- data.frame(tuk2.fam.v$Var1)
tuktab.fam.v <- tuktab.fam.v %>% rownames_to_column("family") 
tuktab.fam.v <- separate(tuktab.fam.v, family, into = c("fam1", "fam2"), sep = "-")
tuktab.fam.v$diff <- NULL
tuktab.fam.v$lwr <- NULL 
tuktab.fam.v$upr<- NULL 
# Make into "wide" format using Var2
tuktab.fam.v <- spread(data = tuktab.fam.v, key = 'fam2', value = 'p.adj', fill = '')
tuktab.fam.v

#And last but not least, Order visuals only
taxa.ord.v <- (data.frame(table(site_visuals$order, site_visuals$date, site_visuals$walk)))
qqnorm(taxa.ord.v$Freq)

shapiro.test(taxa.ord.v$Freq)

## 
##  Shapiro-Wilk normality test
## 
## data:  taxa.ord.v$Freq
## W = 0.10398, p-value < 2.2e-16

#it's normal

ordav <- aov(Freq ~ Var1, taxa.ord.v)
summary(ordav)

##               Df Sum Sq Mean Sq F value Pr(>F)    
## Var1           7    506   72.22   22.07 <2e-16 ***
## Residuals   2072   6779    3.27                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#significant model

#Now see what stood out and in what direction
#Tukey HSD
tuk2.ord.v <- TukeyHSD(ordav, "Var1", ordered = TRUE)
#Make viewable
tuktab.ord.v <- data.frame(tuk2.ord.v$Var1)
tuktab.ord.v <- tuktab.ord.v %>% rownames_to_column("order") 
tuktab.ord.v <- separate(tuktab.ord.v, order, into = c("order1", "order2"), sep = "-")
tuktab.ord.v$diff <- NULL
tuktab.ord.v$lwr <- NULL 
tuktab.ord.v$upr<- NULL 
# Make into "wide" format using Var2
tuktab.ord.v <- spread(data = tuktab.ord.v, key = 'order2', value = 'p.adj', fill = '')
tuktab.ord.v

#Does mesohabitat location differ significantly?
meso.local <- (data.frame(table(site$mesohabitat, site$date, site$walk)))
qqnorm(meso.local$Freq)

shapiro.test(meso.local$Freq)

## 
##  Shapiro-Wilk normality test
## 
## data:  meso.local$Freq
## W = 0.16689, p-value < 2.2e-16

#it's normal

mesoa <- aov(Freq ~ Var1, meso.local)
summary(mesoa)

##               Df Sum Sq Mean Sq F value Pr(>F)    
## Var1           9    461   51.20   14.17 <2e-16 ***
## Residuals   2590   9361    3.61                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#significant model

#Now see what stood out and in what direction
#Tukey HSD
tuk3.meso <- TukeyHSD(mesoa, "Var1", ordered = TRUE)
#Make viewable
tuktab.meso <- data.frame(tuk3.meso$Var1)
tuktab.meso <- tuktab.meso %>% rownames_to_column("mesohabitat") 
tuktab.meso <- separate(tuktab.meso, mesohabitat, into = c("meso1", "meso2"), sep = "-")
tuktab.meso$diff <- NULL
tuktab.meso$lwr <- NULL 
tuktab.meso$upr<- NULL 
# Make into "wide" format using Var2
tuktab.meso <- spread(data = tuktab.meso, key = 'meso2', value = 'p.adj', fill = '')
tuktab.meso

#Compare these results to bargraphs, and the differences are quite apparent! Hard to work with models that have this many variables, but can be combed through. 

Bird Visitation

In this section, will explore how birds visited the cacti with different floral output and size, as well as bird diversity and behavior at cacti.

#bring in datasets of interest
visit <- read.csv("~/Masters/Flower_CactusBirdInteractionR/data/mimic_flower_models/focal_bird_visitation.csv")
head(visit)

cam.trap.specs <- read.csv("~/Masters/Flower_CactusBirdInteractionR/data/mimic_flower_models/cam_trap_specs.csv")

#Visitations (including non-pollinating behaviors)
avg_treatment <- count(visit, vars=treatment, wt_var=experiment)
#fifiteen on flowers experiment did not have any sort of visitation at all
avg_treatment

#basemap for sunset cove
cali <- get_map(location = c(lon = -115.663, lat = 34.7825), zoom = 18, color="bw")

## Source : https://maps.googleapis.com/maps/api/staticmap?center=34.7825,-115.663&zoom=18&size=640x640&scale=2&maptype=terrain&language=en-EN&key=xxx-6x1tmeuVFVYhc8MMkmuUuKOALLNc

#How does mimic model compare to real-flower model in terms of visitation?
ggplot(visit, aes(x=experiment)) + geom_bar() 

#Let's map the location of all cacti in experiments
experiments <- ggmap(cali)
experiments <- experiments + 
  geom_point(data = visit, aes(x = longitude, y = latitude, colour = experiment, group = experiment), alpha = 3/10, size = 4) + 
  labs(title = "Experimental Cacti", x = "longitude", y = "latitude", color = "Experiment")
experiments

#remember, we began using the same cacti for the second half of the experiment.


#How does visitation differ between treatments by experiment type?
ggplot(visit, aes(x=treatment)) + geom_bar() + facet_grid(. ~experiment)

#Notice how pole was the same between both experiments

#Treatment distribution by experiment
treatments <- ggmap(cali)
treatments <- treatments +
  geom_point(data=visit, aes(x=longitude, y=latitude, colour = treatment, group = treatment), alpha = 3/10, size =4) +
  labs(title = "Cacti Location for Each Experiment", x = "longitude", y = "latitude", color = "Treatment") + facet_grid(.~experiment)
treatments

#problem because we don't actually see real-flower treatment of 15 because it had no visitations at all, so no records within this dataset. 


#Distribution of bird species by experiment
ggplot(visit, aes(x=taxa)) + geom_bar() + facet_grid(experiment ~.) + theme(axis.text.x = element_text(angle = 60, hjust =1))

#Map species of bird visitations to the cacti.
birds.at.cacti <- ggmap(cali)
birds.at.cacti <- birds.at.cacti +
  geom_point(data=visit, aes(x=longitude, y=latitude, colour = taxa, group = taxa), alpha = 3/10, size =4) +
  labs(title = "Bird Visitations to Cacti", x = "longitude", y = "latitude", color = "Species") + facet_grid(.~experiment)
birds.at.cacti

#Distribution of behaviors 
ggplot(visit, aes(x=behavior)) + geom_bar() + facet_grid(experiment ~.) + theme(axis.text.x = element_text(angle = 60, hjust =1))

#Interesting to think about behaviors that aren't necessarily pollinating, but are indicative of a bird occupying that space/claiming territory (sparring, calling for example).

#Behavior map
behavior.map <- ggmap(cali)
behavior.map <- behavior.map +
  geom_point(data=visit, aes(x=longitude, y=latitude, colour = behavior, group = behavior), alpha = 3/10, size =4) +
  labs(title = "Cacti Location for Each Experiment", x = "longitude", y = "latitude", color = "Behavior") + facet_grid(.~experiment)
behavior.map

#All we've mapped so far is including non-pollinating behaviors.
#Let's double check for only pollinators
pollinators <- filter(visit, visit$behavior == "Pollinating")
#Notice it removes half of observations! Goes from 88 to 44 observations

ggplot(pollinators, aes(x=treatment)) + geom_bar() + facet_grid(. ~experiment)

#removes treatment zero (no occurrences in either experiment of pollinating at a plant without flowers, this much is obviously expected), and shows no difference between fifteen, pole, and thirty for mimics.

#Pollinating behaviors by treatment for each experiment
pollinate.treatments <- ggmap(cali)
pollinate.treatments <- pollinate.treatments +
  geom_point(data=visit, aes(x=longitude, y=latitude, colour = treatment, group = treatment), alpha = 3/10, size =4) +
  labs(title = "Pollination by treatment", x = "longitude", y = "latitude", color = "Treatment") + facet_grid(.~experiment)
pollinate.treatments

#Map specs for camera traps, which does not take into account redos
#Experiment locations together
cam.traps.ex <- ggmap(cali)
cam.traps.ex <- cam.traps.ex +
  geom_point(data = cam.trap.specs, aes(x=longitude, y=latitude, colour = experiment, group = experiment), alpha = 3/10, size = 4) +
  labs(title = "Cactus Individuals by Experiment", x = "longitude", y = "latitude", color = "Experiment") 
cam.traps.ex

#Treatment locations facetted by experiment
cam.traps.treat <- ggmap(cali)
cam.traps.treat <- cam.traps.treat +
  geom_point(data = cam.trap.specs, aes(x=longitude, y=latitude, colour = treatment, group = treatment), alpha =3/10, size = 4) +
  labs(title = "Treatment Locations: Camera traps", x = "longitude", y = "latitude", color = "Treatment") +
  facet_grid(.~experiment)
cam.traps.treat

#Bring in relevant datasets
hobo <- read.csv("~/Masters/Flower_CactusBirdInteractionR/data/mimic_flower_models/cactus_architecture_temperature.csv")

hobo.specs <- read.csv("~/Masters/Flower_CactusBirdInteractionR/data/mimic_flower_models/pendant_specs.csv")
head(hobo.specs)

#Join hobo.specs to hobo file
hobo <- left_join(hobo, hobo.specs, by = "cactus")
hobo$X <- NULL

#remove na (from cactus 6 and 8)
hobo <- na.omit(hobo, col = c(5:6))

#Map the cacti
arch.map <- ggmap(cali)
arch.map <- arch.map +
  geom_point(data=hobo.specs, aes(x=lon, y=lat, colour = experiment, group = experiment), alpha = 3/10, size =4) +
  labs(title = "Cactus architecture influence on temperature/light", x = "longitude", y = "latitude", color = "experiment") 
arch.map

#remove na (from cactus 6 and 8)
hobo <- na.omit(hobo, col = c(5:6))

#Temperature throughout both experiments, raw data
temp <- ggplot(hobo, aes(x = date.time, y = temp, colour = experiment)) +
  geom_line() +
  theme(axis.text.x = element_blank()) +
  labs(title = "Temperature through both experiments", x = "Time", y = "Tempearture (F)", color = "Experiment")
temp 

#Light throughout both experiments, raw dat
light <- ggplot(hobo, aes(x = date.time, y = light, colour = experiment)) +
  geom_line() +
  theme(axis.text.x = element_blank()) +
  labs(title = "Light intensity through both experiments", x = "Time", y = "Light (Lumens/ft^2)", color = "Experiment")
light

#Get stat summaries by cactus in neat dataframe
summary1 <- hobo %>% group_by(cactus) %>% summarise_at(vars(temp), funs(mean, sd, min, max))
summary2 <- hobo %>% group_by(cactus) %>% summarise_at(vars(light), funs(mean, sd, min, max))
summary <- left_join(summary1, summary2, by = "cactus")
names(summary) <- c("cactus", "mean.temp", "sd.temp", "min.temp", "max.temp", "mean.light", "sd.light", "min.light", "max.light")
summary <- mutate(summary, sd.min.temp = mean.temp - sd.temp) 
summary <-  mutate(summary, sd.max.temp = mean.temp + sd.temp)
summary <-  mutate(summary, sd.min.light = mean.light - sd.light)
summary <-  mutate(summary, sd.max.light = mean.light + sd.light)

#mean temperatures of each cactus
#first as a bargraph
temp.mean.bar <- ggplot(summary, aes(x = cactus, y = mean.temp)) +
  geom_bar(stat = "identity") + geom_errorbar(summary, mapping = aes(ymin = sd.min.temp, ymax = sd.max.temp)) +
  labs(title = "Mean temperatures of each cactus", x = "Cactus ID", y = "Mean Temperature (F)") +
  scale_x_continuous(breaks=seq(1, 20, 1)) +
  theme_minimal()
temp.mean.bar

#then as a boxplot
temp.mean.boxplot <- ggplot(hobo, aes(x = cactus, y = temp, group = cactus)) +
  geom_boxplot() +
  labs(title = "Mean temperatures of each cactus", x = "Cactus ID", y = "Mean temperature (F)") +
  theme_minimal()
temp.mean.boxplot

#mean light intensity for each cactus
light.mean.bar <- ggplot(summary, aes(x = cactus, y = mean.light)) +
  geom_bar(stat = "identity") + geom_errorbar(summary, mapping = aes(ymin = sd.min.temp, ymax = sd.max.temp)) + 
  labs(title = "Mean light intensity of each cactus", x = "Cactus ID", y = "Mean Light Intensity (Lumns/ft^2)") + 
  scale_x_continuous(breaks=seq(1, 20, 1)) +
  theme_minimal()
light.mean.bar

#then as a boxplot
light.mean.boxplot <- ggplot(hobo, aes(x = cactus, y = light, group = cactus)) +
  geom_boxplot() +
  labs(title = "Mean light intensity of each cactus", x = "Cactus ID", y = "Mean temperature (F)") +
  theme_minimal()
light.mean.boxplot

#Get stat summaries by architecture in neat dataframe
summary3 <- hobo %>% group_by(architecture) %>% summarise_at(vars(temp), funs(mean, sd, min, max))
summary4 <- hobo %>% group_by(architecture) %>% summarise_at(vars(light), funs(mean, sd, min, max))
summary.arch <- left_join(summary3, summary4, by = "architecture")
names(summary.arch) <- c("architecture", "mean.temp", "sd.temp", "min.temp", "max.temp", "mean.light", "sd.light", "min.light", "max.light")
summary.arch <- mutate(summary.arch, sd.min.temp = mean.temp - sd.temp) 
summary.arch <-  mutate(summary.arch, sd.max.temp = mean.temp + sd.temp)
summary.arch <-  mutate(summary.arch, sd.min.light = mean.light - sd.light)
summary.arch <-  mutate(summary.arch, sd.max.light = mean.light + sd.light)

#mean temperatures of each architecture type
#first as a bargraph
arch.temp.mean.bar <- ggplot(summary.arch, aes(x = architecture, y = mean.temp)) +
  geom_bar(stat = "identity") + geom_errorbar(summary.arch, mapping = aes(ymin = sd.min.temp, ymax = sd.max.temp)) +
  labs(title = "Mean temperatures of each architecture type", x = "Architecture Type", y = "Mean Temperature (F)") +
  theme_minimal()
arch.temp.mean.bar

#then as a boxplot
arch.temp.mean.boxplot <- ggplot(hobo, aes(x = architecture, y = temp, group = architecture)) +
  geom_boxplot() +
  labs(title = "Mean temperatures of each architecture type", x = "Architecture Type", y = "Mean temperature (F)") +
  theme_minimal()
arch.temp.mean.boxplot

#mean light intensity for each architecture type
arch.light.mean.bar <- ggplot(summary.arch, aes(x = architecture, y = mean.light)) +
  geom_bar(stat = "identity") + geom_errorbar(summary.arch, mapping = aes(ymin = sd.min.temp, ymax = sd.max.temp)) + 
  labs(title = "Mean light intensity of each architecture type", x = "Architecture Type", y = "Mean Light Intensity (Lumns/ft^2)") + 
  theme_minimal()
arch.light.mean.bar

#doesn't look great because the standard deviation is so much lower (there are some extremely high "outliers", in the statistical sense. This is to be expected, since it's the desert.)

#then as a boxplot
arch.light.mean.boxplot <- ggplot(hobo, aes(x = architecture, y = light, group = architecture)) +
  geom_boxplot() +
  labs(title = "Mean light intensity of each architecture type", x = "Architecture Type", y = "Mean temperature (F)") +
  theme_minimal()
arch.light.mean.boxplot

#Some simple modelling

#Is it normal? 
#Can't use shaprio test because it's over 5000 samples, so let's just do some normal q-q plots
qqnorm(hobo$temp)

#looks pretty normal
qqnorm(hobo$light)

#looks pretty normal
#Good to go ahead with parametric tests

#test for temperature and light correlation
cor.test(hobo$temp, hobo$light, method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  hobo$temp and hobo$light
## t = 206.07, df = 101050, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5396135 0.5482960
## sample estimates:
##       cor 
## 0.5439693

#They are indeed correlated (p < 2.2e-16 & R = 0.54, so they are positively correlated)

#Visualize correlation
temp.light.cor <- ggscatter(hobo, x = "temp", y = "light", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Temperature (F)", ylab = "Light Intensity (lumens/ft^2")
temp.light.cor

#yikes not very pretty, let's try it with means
mean.cor <- ggscatter(summary, x = "mean.temp", y = "mean.light", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Temperature (F)", ylab = "Light Intensity (lumens/ft^2")
mean.cor

#Interestingly, the means are not correlated. Actually, that makes perfect sense...

#Test for difference between architecture types
#Temperature
arch.temp.anova <- aov(hobo$temp ~ hobo$architecture)
summary(arch.temp.anova)

##                       Df   Sum Sq Mean Sq F value Pr(>F)    
## hobo$architecture      2   659275  329638   860.6 <2e-16 ***
## Residuals         101048 38705839     383                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#big significance! p < 2e-16. Interesting, because the boxplot made them look not very different
arch.temp.tukey <- TukeyHSD(arch.temp.anova)
arch.temp.tukey

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = hobo$temp ~ hobo$architecture)
## 
## $`hobo$architecture`
##                     diff       lwr       upr p adj
## middle-bottom -5.7658777 -6.119345 -5.412411 0e+00
## top-bottom    -4.9877047 -5.341177 -4.634232 0e+00
## top-middle     0.7781729  0.424753  1.131593 7e-07

#Light Intensity
arch.light.anova <- aov(hobo$light ~ hobo$architecture)
summary(arch.light.anova)

##                       Df    Sum Sq  Mean Sq F value Pr(>F)    
## hobo$architecture      2 1.382e+10 6.91e+09   511.7 <2e-16 ***
## Residuals         101048 1.365e+12 1.35e+07                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Again, big significance! p < 2e-16. 
arch.light.tukey <- TukeyHSD(arch.light.anova)
arch.light.tukey

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = hobo$light ~ hobo$architecture)
## 
## $`hobo$architecture`
##                    diff       lwr       upr   p adj
## middle-bottom -135.9670 -202.3353 -69.59861 4.7e-06
## top-bottom     707.5701  641.2008 773.93946 0.0e+00
## top-middle     843.5371  777.1776 909.89658 0.0e+00

#Test for influence of cactus size (individual metrics and volume as a whole) on temperature for each architecture type
#ANCOVA because explanatory variables are both categorical (architecture type) and continuous (size)
arch.volume.temp.model <- lm(temp ~ architecture * volume, data = hobo)
summary(arch.volume.temp.model)

## 
## Call:
## lm(formula = temp ~ architecture * volume, data = hobo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -96.698 -14.987  -3.044  14.004  84.049 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                7.621e+01  1.351e-01 563.888   <2e-16 ***
## architecturemiddle        -5.762e+00  1.911e-01 -30.152   <2e-16 ***
## architecturetop           -5.255e+00  1.911e-01 -27.497   <2e-16 ***
## volume                    -5.346e-07  2.622e-08 -20.389   <2e-16 ***
## architecturemiddle:volume -1.267e-09  3.708e-08  -0.034   0.9727    
## architecturetop:volume     8.347e-08  3.708e-08   2.251   0.0244 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.46 on 101045 degrees of freedom
## Multiple R-squared:  0.02762,    Adjusted R-squared:  0.02757 
## F-statistic: 573.9 on 5 and 101045 DF,  p-value: < 2.2e-16

#ISSUE? Should "architecture" be split into wide format and each variable (top, middle and bottom) be treated as an indivdiual covariate?
#Intercept here represents architecturebottom, because it was left out as the reference variable?
#P-values tell us models were significant. Looking at the coefficients estimates, it appears that the middle is about 6 degree cooler than the bottom, and the top is about 5 degree cooler. But when you consider volume of the cactus as an influencing factor (:) then the middle is only 1 degree cooler than the bottom (but this model was not significant so ignore it), and the top is 8 degrees hotter (This makes so much sense!) Bigger cacti provide more shade to the bottom, smaller plants don't!!

Cactus allocation

#bring in dataset
allo <- read.csv("~/Masters/Flower_CactusBirdInteractionR/data/paired_flower_fruit/paired_flower_fruit.csv")
head(allo)

#basemap closer to survey site
paired.fruit.flower.basemap <- get_map(location = c(lon = -115.662, lat = 34.7813), zoom = 19, color = "bw")

## Source : https://maps.googleapis.com/maps/api/staticmap?center=34.7813,-115.662&zoom=19&size=640x640&scale=2&maptype=terrain&language=en-EN&key=xxx-6x1tmeuVFVYhc8MMkmuUuKOALLNc

#Bubble map of total plant volume compared to buds
allo.volume.buds <- ggmap(paired.fruit.flower.basemap)
allo.volume.buds <- allo.volume.buds +
   scale_color_gradient(low="#8CF554", high="#F59154") +
   geom_point(data = allo, aes(x = lon, y = lat, size = volume, color = buds), alpha = 6/10) 
allo.volume.buds

#Bubble map of new growth compared to buds
allo.growth.buds <- ggmap(paired.fruit.flower.basemap)
allo.growth.buds <- allo.growth.buds + 
  scale_color_gradient(low = "#8CF554", high="#F59154") +
  geom_point(data = allo, aes(x = lon, y = lat, size = newgrowth, color = buds), alpha = 6/10) 
allo.growth.buds

#Some means to wrap our brains around the numbers
mean(allo$volume)

## [1] 147117.9

mean(allo$buds)

## [1] 120.89

mean(allo$newgrowth)

## [1] 27.71

#Normality?
shapiro.test(allo$volume) #normal

## 
##  Shapiro-Wilk normality test
## 
## data:  allo$volume
## W = 0.72487, p-value = 2.132e-12

shapiro.test(allo$newgrowth) #normal

## 
##  Shapiro-Wilk normality test
## 
## data:  allo$newgrowth
## W = 0.7429, p-value = 6.103e-12

shapiro.test(allo$buds) #normal

## 
##  Shapiro-Wilk normality test
## 
## data:  allo$buds
## W = 0.76037, p-value = 1.776e-11

#All normal, moving on...

#See some correlations
allo.cor.plot1 <- ggplot(data = allo, aes(x=volume, y=newgrowth)) + geom_point() + geom_smooth()
allo.cor.plot1

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

allo.cor1 <- cor.test(allo$volume, allo$newgrowth, method = 'pearson')
allo.cor1

## 
##  Pearson's product-moment correlation
## 
## data:  allo$volume and allo$newgrowth
## t = 5.1735, df = 98, p-value = 1.22e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2934492 0.6045865
## sample estimates:
##       cor 
## 0.4631707

#Yes, positively correlated (R = .46, p= 1.22e-06)

allo.cor.plot2 <- ggplot(data = allo, aes(x=volume, y=buds)) + geom_point() + geom_smooth()
allo.cor.plot2

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

allo.cor2 <- cor.test(allo$volume, allo$buds, method = 'pearson')
allo.cor2

## 
##  Pearson's product-moment correlation
## 
## data:  allo$volume and allo$buds
## t = 6.4685, df = 98, p-value = 3.904e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3927785 0.6712911
## sample estimates:
##       cor 
## 0.5469965

#Yes, positively correlated (R = .55, p = 3.904e-09)

allo.cor.plot3 <- ggplot(data = allo, aes(x=buds, y=newgrowth)) + geom_point() + geom_smooth()
allo.cor.plot3

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

allo.cor3 <- cor.test(allo$newgrowth, allo$buds, method = 'pearson')
allo.cor3

## 
##  Pearson's product-moment correlation
## 
## data:  allo$newgrowth and allo$buds
## t = 4.3997, df = 98, p-value = 2.76e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2278939 0.5580352
## sample estimates:
##       cor 
## 0.4061325

#Hmm, they are also positively correlated (R = .40, p = 2.76e-5). I think we need to consider that a bigger plant has more branches, so unless you consider the interaction between newgrowth (or buds) and volume, you won't have something that shows us what's really happening.

#Model for how new growth/volume impacts bud growth
allo.m1 <- lm(buds ~ newgrowth:volume, data=allo)
summary(allo.m1)

## 
## Call:
## lm(formula = buds ~ newgrowth:volume, data = allo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -192.53  -91.13  -54.71   20.13  501.63 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      9.112e+01  1.551e+01   5.876 5.79e-08 ***
## newgrowth:volume 4.677e-06  9.734e-07   4.805 5.57e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 142.2 on 98 degrees of freedom
## Multiple R-squared:  0.1907, Adjusted R-squared:  0.1824 
## F-statistic: 23.09 on 1 and 98 DF,  p-value: 5.574e-06

#Model1 is significant, and the interaction between newgrowth by volume is lowered by about 5 for each addition of a bud
allo.m2 <- lm(buds ~ newgrowth*volume, data = allo)
summary(allo.m2)

## 
## Call:
## lm(formula = buds ~ newgrowth * volume, data = allo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -241.08  -60.89  -30.91    4.52  424.01 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       4.398e+00  2.482e+01   0.177  0.85976    
## newgrowth         2.216e+00  8.134e-01   2.725  0.00764 ** 
## volume            5.308e-04  1.133e-04   4.685 9.24e-06 ***
## newgrowth:volume -3.617e-06  2.033e-06  -1.779  0.07842 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 128.7 on 96 degrees of freedom
## Multiple R-squared:  0.3503, Adjusted R-squared:   0.33 
## F-statistic: 17.26 on 3 and 96 DF,  p-value: 4.85e-09

#This model may be overfit, as newgrowth and volume individually will obviously not be linked with bud production, they have to be together.

#Model for how buds/volume impacts new growth
allo.m3 <- lm(newgrowth ~ buds:volume, data = allo)
summary(allo.m3)

## 
## Call:
## lm(formula = newgrowth ~ buds:volume, data = allo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -57.453 -12.635  -7.211   7.369  87.506 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.221e+01  2.671e+00   8.316 5.34e-13 ***
## buds:volume 1.633e-07  3.381e-08   4.830 5.03e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24.16 on 98 degrees of freedom
## Multiple R-squared:  0.1923, Adjusted R-squared:  0.1841 
## F-statistic: 23.33 on 1 and 98 DF,  p-value: 5.034e-06

#Model is significant, and the itneraction between buds by volume raises the number of new growth by 1.633e07 (so almost 0).
allo.m4 <- lm(newgrowth ~ buds*volume, data = allo)
summary(allo.m4)

## 
## Call:
## lm(formula = newgrowth ~ buds * volume, data = allo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -46.817 -13.444  -3.466   7.273  83.234 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.535e+01  3.655e+00   4.199    6e-05 ***
## buds         4.194e-02  2.296e-02   1.827   0.0709 .  
## volume       5.564e-05  2.406e-05   2.313   0.0229 *  
## buds:volume -2.646e-08  7.797e-08  -0.339   0.7351    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.54 on 96 degrees of freedom
## Multiple R-squared:  0.2487, Adjusted R-squared:  0.2253 
## F-statistic: 10.59 on 3 and 96 DF,  p-value: 4.418e-06

#Now, the model is not significant, but the coefficient is negative now, but still by almost 0. 

#So it seems like the most sensical direction is buds being the dependant variable

Some takeaways to consider for August:

Should we continue camera trapping, focal observations, and cube video-ing? I’m inclined to say scrap focal observations and instead up the frequency of video footage and camera trapping. Although camera traps are not picking up hummingbirds accurately, more comprehensive data on other, non-bird visitors would still be a nice dataset. Camera trapping is also very low effort when raised above grass-line.

Should we deploy cameras/cubes at shrub and open areas still? I’m inclined to say no. Originally, these were options as alternate controls in my design, and we scrapped them. I brought them back because they were relatively low effort, but I still do not think these were necessary or answering any precise questions that the other four treatments did not address.Perhaps best to scrap again for August and save time.